Entrega I RESUELTA

Author

C. Tangana - DNI: 00000000-X

Instrucciones (leer antes de empezar)

  • Modifica de la cabecera del documento .qmd tus datos personales (nombre y DNI). IMPORTANTE: no modifiques nada más de la cabecera.

  • Asegúrate, ANTES de seguir editando el documento, que el archivo .qmd se renderiza correctamente y se genera el .html correspondiente en tu carpeta local de tu ordenador (pulsa el botón Render o guarda el documento con Render on Save activado).

  • Los chunks (cajas de código) creados están o vacíos o incompletos, de ahí que la mayoría tengan la opción #| eval: false. Una vez que edites lo que consideres, debes ir cambiando cada chunck a #| eval: true (o quitarlo directamente) para que se ejecuten.

  • Recuerda que puedes ejecutar chunk a chunk con el botón play o ejecutar todos los chunk hasta uno dado (con el botón a la izquierda del anterior).

  • IMPORTANTE: solo se podrá subir un archivo .html al campus, no se evaluarán entregas sin el .html correctamente generado. Recuerda: el código es una herramienta, no el fin en esta asignatura. Se evaluará especialmente las interpretaciones y conclusiones que detalles. Un ejercicio con el código perfecto pero sin ningún tipo de razonamiento, interpretación o conclusión no superará el 4 sobre 10 de nota.

Paquetes necesarios

Necesitaremos los siguientes paquetes (haz play en el chunk para que se carguen):

rm(list = ls()) # Borramos variables de environment

# descomentar si es la primera vez (y requieren instalación)
# install.packages("tidyverse")
# install.packages("performance")
# install.packages("olsrr")
# install.packages("corrr")
# install.packages("corrplot")
# install.packages("skimr")
library(tidyverse)
library(performance)
library(olsrr)
library(corrr)
library(corrplot)
library(skimr)

Caso práctico I: análisis de datos de cáncer

El archivo de datos a usar lo cargaremos desde el csv cancer.csv.

# no cambies código
datos <- read_csv(file = "./cancer.csv")
datos
# A tibble: 3,047 × 13
   deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge Geography
       <dbl>     <dbl>      <dbl>          <dbl>       <dbl>     <dbl> <chr>    
 1      165.     61898     260131           11.2       500.       39.3 Kitsap C…
 2      161.     48127      43269           18.6        23.1      33   Kittitas…
 3      175.     49348      21026           14.6        47.6      45   Klickita…
 4      195.     44243      75882           17.1       343.       42.8 Lewis Co…
 5      144.     49955      10321           12.5         0        48.3 Lincoln …
 6      176      52313      61023           15.6       180.       45.4 Mason Co…
 7      176.     37782      41516           23.2         0        42.6 Okanogan…
 8      184.     40189      20848           17.8         0        51.7 Pacific …
 9      190.     42579      13088           22.3         0        49.3 Pend Ore…
10      178.     60397     843954           13.1       428.       35.8 Pierce C…
# ℹ 3,037 more rows
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
#   PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
#   BirthRate <dbl>

Los datos representan diferentes características socioeconómicas de distintas regiones, extraídas de the American Community Survey (census.gov), clinicaltrials.gov, y cancer.gov.

¿Nuestro objetivo? Ser capaces de predecir la variable deathRate, nuestra variable objetivo, que representa la mortalidad media de cancer, por cada 100 000 habitantes. La idea es que seamos capaz de determinar cual de las variables es la que más efecto (lineal) tiene en la mortalidad de cáncer, realizando el ajuste posterior.

El resto de variables son:

  • medianIncome: mediana de los ingresos de la región.
  • popEst2015: población de la región
  • povertyPercent: porcentaje de población en situación de pobreza.
  • studyPerCap: ensayos clínicos relacionados por el cáncer realizados por cada 100 000 habitantes.
  • MedianAge: edad mediana.
  • Geography: nombre de la región.
  • AvgHouseholdSize: tamaño medio de los hogares.
  • PercentMarried: porcentaje de habitantes casados.
  • PctHS25_Over: porcentaje de residentes por encima de los 25 años con (máximo) título de bachillerato.
  • PctUnemployed16_Over: porcentaje de residentes de más de 16 años en paro.
  • PctPrivateCoverage: porcentaje de residentes con seguro de salud privado.
  • BirthRate: tasa de natalidad.

IMPORTANTE: siempre que puedas, relaciona los valores numéricos de la salida de R con su fórmula.

Ejercicio 1 (1/10)

Ejecuta el código que consideres para chequear que a) todas las variables son numéricas; b) no hay datos ausentes; c) no hay problemas de rango/codificación. En caso de que alguna de ellas no se cumpla, ejecuta el código que consideres para corregirlo y/o eliminar las variables que consideres. Sé conciso en el código.

Simplemente había que darse cuenta de que

  1. No había ausentes (no se hace nada)
  2. Geography era no numérica –> fuera (de momento).
  3. los rangos permitidos:
  • deathRate: entre 0 e infinito.
  • medianIncome: entre 0 e infinito.
  • popEst2015: entre 0 y 8.000.000.000 (habitantes del mundo xd)
  • povertyPercent: entre 0 y 100 (al ser porcentaje)
  • studyPerCap: entre 0 e infinito
  • MedianAge: entre 0 y 130 años (por poner un tope)
  • AvgHouseholdSize: entre 0 e infinito (si existe alguna finca enorme)
  • PercentMarried: entre 0 y 100 (al ser porcentaje)
  • PctHS25_Over: entre 0 y 100 (al ser porcentaje)
  • PctUnemployed16_Over: entre 0 y 100 (al ser porcentaje)
  • ctPrivateCoverage: entre 0 y 100 (al ser porcentaje)
  • BirthRate: entre 0 e infinito.

Lo comprobamos con skim() o summary()

datos |> skim()
Data summary
Name datos
Number of rows 3047
Number of columns 13
_______________________
Column type frequency:
character 1
numeric 12
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
Geography 0 1 16 42 0 3047 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
deathRate 0 1 178.66 27.75 59.70 161.20 178.10 195.20 362.80 ▁▇▆▁▁
medIncome 0 1 47063.28 12040.09 22640.00 38882.50 45207.00 52492.00 125635.00 ▇▇▁▁▁
popEst2015 0 1 102637.37 329059.22 827.00 11684.00 26643.00 68671.00 10170292.00 ▇▁▁▁▁
povertyPercent 0 1 16.88 6.41 3.20 12.15 15.90 20.40 47.40 ▃▇▃▁▁
studyPerCap 0 1 155.40 529.63 0.00 0.00 0.00 83.65 9762.31 ▇▁▁▁▁
MedianAge 0 1 45.27 45.30 22.30 37.70 41.00 44.00 624.00 ▇▁▁▁▁
AvgHouseholdSize 0 1 2.48 0.43 0.02 2.37 2.50 2.63 3.97 ▁▁▃▇▁
PercentMarried 0 1 51.77 6.90 23.10 47.75 52.40 56.40 72.50 ▁▂▇▇▁
PctHS25_Over 0 1 34.80 7.03 7.50 30.40 35.30 39.65 54.80 ▁▂▇▇▁
PctUnemployed16_Over 0 1 7.85 3.45 0.40 5.50 7.60 9.70 29.40 ▅▇▁▁▁
PctPrivateCoverage 0 1 64.35 10.65 22.30 57.20 65.10 72.10 92.30 ▁▂▇▇▂
BirthRate 0 1 5.64 1.99 0.00 4.52 5.38 6.49 21.33 ▂▇▁▁▁

La única con problemas es MedianAge que tiene de máximo 600 años: filtramos por ese umbral de 130 años y los eliminamos (más adelante los pasaremos a NA y los depuraremos).

datos_preproc <-
  datos |>
  select(where(is.numeric)) |> 
  filter(MedianAge < 130)
datos_preproc
# A tibble: 3,017 × 12
   deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
       <dbl>     <dbl>      <dbl>          <dbl>       <dbl>     <dbl>
 1      165.     61898     260131           11.2       500.       39.3
 2      161.     48127      43269           18.6        23.1      33  
 3      175.     49348      21026           14.6        47.6      45  
 4      195.     44243      75882           17.1       343.       42.8
 5      144.     49955      10321           12.5         0        48.3
 6      176      52313      61023           15.6       180.       45.4
 7      176.     37782      41516           23.2         0        42.6
 8      184.     40189      20848           17.8         0        51.7
 9      190.     42579      13088           22.3         0        49.3
10      178.     60397     843954           13.1       428.       35.8
# ℹ 3,007 more rows
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
#   PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
#   BirthRate <dbl>

Ejercicio 2 (0.8/10)

Realiza un primer análisis exploratorio para entender nuestros datos (cada variable por separado). Hazlo tanto de manera numérica como de manera visual (no hace falta que ejecutas 19 000 gráficos, decide uno de los aprendidos para ello y saca conclusiones, al grano). ¿Hay datos atípicos?

Por ejemplo densidades…

datos_tidy <-
  datos_preproc |> 
  pivot_longer(cols = everything(), names_to = "variable", values_to = "values")

ggplot(datos_tidy) + 
  geom_density(aes(x = values, color = variable, fill = variable)) +
  facet_wrap(~variable, scale = "free") +
  theme_minimal()

…o boxplots

ggplot(datos_tidy) + 
  geom_boxplot(aes(x = values, color = variable, fill = variable)) +
  facet_wrap(~variable, scale = "free") +
  theme_minimal()

Si se observan atípicos, en especial en popEst2015 (hay algunas regiones muy grandes), en studyPerCap (es decir, no en todas las regiones se investiga por igual, hay algunas donde se hacen muchos estudios), y las variables de renta (medIncome y PovertyPercent), que muestran una gran desigualdad.

De momento suficiente con darse cuenta de esto, ya intervendremos.

Ejercicio 3 (1.2/10)

Nuestro objetivo será predecir la mortalidad de cáncer. De momento trabajaremos en el contexto de la regresión univariante, así que lo primero que haremos tras preprocesar y entender los datos será seleccionar la variable predictora que mayor efecto lineal tenga respecto a la variable objetivo. Hazlo tanto de manera numérica como visual, de manera que puedas descartar correlaciones espúreas (que no haya un dinosaurio, por ejemplo).

Sacamos primero la matriz de correlaciones

datos_preproc |> correlate()
Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'
# A tibble: 12 × 13
   term      deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
   <chr>         <dbl>     <dbl>      <dbl>          <dbl>       <dbl>     <dbl>
 1 deathRate  NA         -0.428     -0.119          0.429     -0.0225   -0.00429
 2 medIncome  -0.428     NA          0.235         -0.788      0.0444   -0.117  
 3 popEst20…  -0.119      0.235     NA             -0.0649     0.0557   -0.176  
 4 povertyP…   0.429     -0.788     -0.0649        NA         -0.0562   -0.194  
 5 studyPer…  -0.0225     0.0444     0.0557        -0.0562    NA        -0.0317 
 6 MedianAge  -0.00429   -0.117     -0.176         -0.194     -0.0317   NA      
 7 AvgHouse…  -0.0381     0.113      0.109          0.0736    -0.00388  -0.357  
 8 PercentM…  -0.267      0.354     -0.161         -0.642     -0.0383    0.430  
 9 PctHS25_…   0.405     -0.471     -0.312          0.194     -0.0852    0.331  
10 PctUnemp…   0.380     -0.452      0.0519         0.654     -0.0314   -0.128  
11 PctPriva…  -0.385      0.724      0.0525        -0.822      0.0934    0.0692 
12 BirthRate  -0.0876    -0.0109    -0.0573        -0.0123     0.0107   -0.101  
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
#   PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
#   BirthRate <dbl>

Con focus() del paquete {corrr} podemos filtrar solo la que nos interesa: predictoras vs objetivo

datos_preproc |>
  correlate() |> 
  corrr::focus(deathRate)
Correlation computed with
• Method: 'pearson'
• Missing treated using: 'pairwise.complete.obs'
# A tibble: 11 × 2
   term                 deathRate
   <chr>                    <dbl>
 1 medIncome             -0.428  
 2 popEst2015            -0.119  
 3 povertyPercent         0.429  
 4 studyPerCap           -0.0225 
 5 MedianAge             -0.00429
 6 AvgHouseholdSize      -0.0381 
 7 PercentMarried        -0.267  
 8 PctHS25_Over           0.405  
 9 PctUnemployed16_Over   0.380  
10 PctPrivateCoverage    -0.385  
11 BirthRate             -0.0876

Y podemos visualizarlas

datos_preproc |>
  cor() |> 
  corrplot(method = "color")

La variable con una mayor correlación es povertyPercent (en positivo, a más pobreza, más mortalidad) y la variable medIncome (en negativo, a más ingresos, menos mortalidad). El signo importa y su interpretación también.

Antes de quedarnos con una de ellas debemos comprobar que no son correlaciones espúreas, así que visualizamos todas vs objetivo

datos_tidy <-
  datos_preproc |> 
  pivot_longer(cols = -deathRate, names_to = "variable", values_to = "values")

ggplot(datos_tidy, aes(x = values, y = deathRate, color = variable)) + 
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~variable, scale = "free") +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

Aunque sea una relación débil, no parece aberrante que podamos ajustar una recta a esa nube de puntos (no hay por ejemplo un…dinosaurio). Nos quedaremos por tanto con povertyPercent como predictora

Ejercicio 4 (1/10)

Una vez elegida la variable con mayor efecto lineal, el objetivo será usar dicha covariable para predecir la variable objetivo. Lo que haremos será separar las observaciones en dos datasets distintos: un dataset train con el que entrenaremos el modelo (el que usará la regresión para aprender y determinar sus coeficientes) y un segundo dataset test que usaremos en el futuro para obtener predicciones. Para ello realiza un muestreo aleatorio de manera que train contenga el 80% de los datos y test el 20% restante. IMPORTANTE: no cambies la semilla.

Primero vamos a crear un id de la fila, para luego poder filtrar correctamente. Para ello usaremos rowid_to_column() que nos numera automáticamente las filas.

Después usaremos slice_sample(prop = 0.8, replace = FALSE). Importante: replace = FALSE viene ya así por defecto, pero es importante que entendamos por qué: no quieroe que haya dos observaciones repetidas.

Y por último, para el otro 20% no puedo hacer otro slice_sample ya que no me garantiza que no caigan observaciones que ya están en train. Así que usaremos un anti_join o un filter()

set.seed(12345)

# Creamos id
datos_preproc <-
  datos_preproc |>
  rowid_to_column(var = "id")
train <- datos_preproc |> slice_sample(prop = 0.8)

# versión antijoin
test <- datos_preproc |> anti_join(train, by = "id")
test
# A tibble: 604 × 13
      id deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
   <int>     <dbl>     <dbl>      <dbl>          <dbl>       <dbl>     <dbl>
 1     5      144.     49955      10321           12.5         0        48.3
 2    25      237.     37625      23372           22.1         0        41.9
 3    39      207.     42121      68714           13.9       160.       42  
 4    49      210.     35799      27037           20.6      1147.       43  
 5    53      156.     43835     104236           22.5        95.9      30  
 6    58      140.     38204       7229           16.7         0        49.1
 7    64      169.     37468      29126           20.8         0        43.2
 8    85      158.     68430      49762            5.9         0        39  
 9   103      138.     55955       5937            9.2         0        38.8
10   109      163.     40979       3625           12.7         0        47.5
# ℹ 594 more rows
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
#   PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
#   BirthRate <dbl>
# version filter
test <- datos_preproc |> filter(!(id %in% train$id))
test
# A tibble: 604 × 13
      id deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
   <int>     <dbl>     <dbl>      <dbl>          <dbl>       <dbl>     <dbl>
 1     5      144.     49955      10321           12.5         0        48.3
 2    25      237.     37625      23372           22.1         0        41.9
 3    39      207.     42121      68714           13.9       160.       42  
 4    49      210.     35799      27037           20.6      1147.       43  
 5    53      156.     43835     104236           22.5        95.9      30  
 6    58      140.     38204       7229           16.7         0        49.1
 7    64      169.     37468      29126           20.8         0        43.2
 8    85      158.     68430      49762            5.9         0        39  
 9   103      138.     55955       5937            9.2         0        38.8
10   109      163.     40979       3625           12.7         0        47.5
# ℹ 594 more rows
# ℹ 6 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
#   PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
#   BirthRate <dbl>

Ejercicio 5 (1.8/10)

No usaremos el dataset test hasta el final para evaluar las predicciones (con un dataset que el modelo no conoce). Con el dataset train a) ejecuta el ajuste de una regresión lineal univariante; b) interpreta los coeficientes; c) realiza la diagnosis de la manera más completa posible; d) comenta todo lo que sepas sobre la segunda tabla de la salida (la de inferencia de los parámetros). Elimina variables si su efecto no fuese significativo.

Tras separar los datos, entrenamos el modelo con train

ajuste <- lm(data = train, formula = deathRate ~ povertyPercent)
ajuste |> summary()

Call:
lm(formula = deathRate ~ povertyPercent, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-121.490  -14.128    1.459   15.520  170.393 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    147.53777    1.44389  102.18   <2e-16 ***
povertyPercent   1.84645    0.07942   23.25   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 25.2 on 2411 degrees of freedom
Multiple R-squared:  0.1831,    Adjusted R-squared:  0.1828 
F-statistic: 540.5 on 1 and 2411 DF,  p-value: < 2.2e-16

De momento lo único que podemos interpretar es

  • \(\beta_0 = 147.53777\), es decir, la predicción del modelo da un ratio de mortalidad de \(147.53777\) (por cada 100 000 habitantes) cuando no hay un solo habitante pobre en la región.
  • \(\beta_1 = 1.84645\), es decir, por cada punto que suba el porcentaje de pobreza, la mortalidad media por cada 100 000 sube \(1.84645\) puntos.

Así el ajuste formulado sería

\[Mortalidad = 147.53777 + 1.84645 * \%pobreza\]

¿Se cumplen las hipótesis?

check_model(ajuste)

A priori parece que gráficamente se cumplen todas. Vamos a chequear una a una

  1. Errores incorrelados: parece existir independencia entre los errores
check_autocorrelation(ajuste)
OK: Residuals appear to be independent and not autocorrelated (p = 0.228).
  1. Normalidad:
ols_test_normality(ajuste)
Warning in ks.test.default(y, "pnorm", mean(y), sd(y)): ties should not be
present for the Kolmogorov-Smirnov test
-----------------------------------------------
       Test             Statistic       pvalue  
-----------------------------------------------
Shapiro-Wilk              0.9834         0.0000 
Kolmogorov-Smirnov        0.0409          6e-04 
Cramer-von Mises         194.4758        0.0000 
Anderson-Darling          7.7665         0.0000 
-----------------------------------------------

Según el p-valor no hay normalidad…pero vamos a asumir que sí por el siguiente gráfico

ggplot(tibble("residuals" = ajuste$residuals)) +
  stat_qq(aes(sample = residuals)) + 
  stat_qq_line(aes(sample = residuals)) +
  theme_minimal()

Los valores se pegan a la línea, sobre todo en el centro que es lo importante.

Además si pintamos su densidad vs una densidad de una normal (con misma media y desviación)

set.seed(12345)
ggplot(tibble("residuals" = ajuste$residuals,
              "normal" =
                rnorm(n = length(ajuste$residuals), mean = mean(ajuste$residuals),
                      sd = sd(ajuste$residuals)))) +
  geom_density(aes(x = residuals), color = "#822181",
               linewidth = 1.5) +
  geom_density(aes(x = normal), color = "#EFE981",
               linewidth = 1.5) +
  theme_minimal()

Vemos que por lo que puede estar fallando la normalidad es que, por tener poco tamaño muestral E IMPORTANTE NO ESTAR TRATANDO OUTLIERS, la distribución de los errores nos sale ligeramente más apuntada hacía arriba (leptocúrtica, kurtosis positiva).

Esto no había que hacerlo de momento, pero para entender porque los p-valores no concuerdan con lo que ven nuestros ojos, vamos a detectar outliers en la variable predictora

Para eso vamos a usar rápido la función scores() del paquete {outliers}: si ponemos en type = "iqr" (es decir, nos vamos a basar en un boxplot para detectar outliers), nos devuelve cuanto se aleja cada dato de los límites del rango intercuartílico

distancia_x <- 
  train$povertyPercent |>
  outliers::scores(type = "iqr")
distancia_x[1:20]
 [1]  0.00000000  0.00000000  2.60975610  0.01219512  0.19512195  0.00000000
 [7]  0.53658537  0.00000000  1.09756098  0.93902439 -0.37804878  0.00000000
[13] -0.53658537  0.00000000  0.00000000  0.43902439  0.00000000  0.00000000
[19]  0.00000000  0.00000000
distancia_y <- 
  train$deathRate |>
  outliers::scores(type = "iqr")
distancia_y[1:20]
 [1]  0.00000000  0.08430233  1.51453488  0.71802326  0.55232558  0.00000000
 [7]  0.00000000 -0.28779070  0.00000000  0.60755814 -0.08139535  0.00000000
[13] -0.90988372  0.00000000  0.61627907  0.00000000 -1.18313953 -1.70058140
[19]  0.00000000 -0.10465116

Un forma típica y rápida de detectar outliers es cuando dicha distancia supera el 1.5 (en valor absoluto), así que lo que hacemos es eliminar dichas observaciones

train_sin_outliers <-
  train |> 
  filter(abs(distancia_x) < 1.5 & abs(distancia_y) < 1.5)

ajuste_outliers <-
  lm(data = train_sin_outliers, formula = deathRate ~ povertyPercent)
ols_test_normality(ajuste_outliers)
Warning in ks.test.default(y, "pnorm", mean(y), sd(y)): ties should not be
present for the Kolmogorov-Smirnov test
-----------------------------------------------
       Test             Statistic       pvalue  
-----------------------------------------------
Shapiro-Wilk              0.9952         0.0000 
Kolmogorov-Smirnov        0.027          0.0676 
Cramer-von Mises         186.1388        0.0000 
Anderson-Darling          2.941          0.0000 
-----------------------------------------------

Empieza a ser normal… Con un tamaño muestral más grande y un tratamiento más rigurosos de los atípicos lo tendríamos.

  1. Homocedasticidad
check_heteroscedasticity(ajuste)
Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).
check_heteroscedasticity(ajuste_outliers)
Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

De la misma manera no parece que se cumpla a priori la homocedasticidad (por lo mismo) pero tenemos lo que nos interesa: errores en torno a una banda constante de varianza

ggplot(tibble("id" = 1:length(ajuste$residuals),
              "residuals" = ajuste$residuals),
       aes(x = id, y = residuals)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

ggplot(tibble("id" = 1:length(ajuste_outliers$residuals),
              "residuals" = ajuste_outliers$residuals),
       aes(x = id, y = residuals)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

En ambos casos están en torno a una banda, en el segundo aún más ya que se ha eliminado uno de los outliers en los residuales. Si siguiésemos con un tratamiento mejor de los mismos, voilá.

  1. Linealidad:
ajuste_res <-
  lm(data = tibble("residuals" = ajuste$residuals,
                   "fitted" = ajuste$fitted.values),
     formula = residuals ~ fitted)
ajuste_res |> summary()

Call:
lm(formula = residuals ~ fitted, data = tibble(residuals = ajuste$residuals, 
    fitted = ajuste$fitted.values))

Residuals:
     Min       1Q   Median       3Q      Max 
-121.490  -14.128    1.459   15.520  170.393 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.394e-14  7.713e+00       0        1
fitted      -1.395e-16  4.301e-02       0        1

Residual standard error: 25.2 on 2411 degrees of freedom
Multiple R-squared:  2.366e-32, Adjusted R-squared:  -0.0004148 
F-statistic: 5.705e-29 on 1 and 2411 DF,  p-value: 1
ajuste_res <-
  lm(data = tibble("residuals" = ajuste$residuals,
                   "fitted" = ajuste$fitted.values),
     formula = residuals ~ fitted + I(fitted^2))
chisq.test(ajuste$residuals, ajuste$fitted.values,
           simulate.p.value = TRUE, B = 300)

    Pearson's Chi-squared test with simulated p-value (based on 300
    replicates)

data:  ajuste$residuals and ajuste$fitted.values
X-squared = 784225, df = NA, p-value = 0.003322

Respecto a la linealidad ídem: aunque el p-valor no salga como debería, no se aprecia patrón (ni lineal ni cuadrático)

ggplot(tibble("residuals" = ajuste$residuals,
                   "fitted" = ajuste$fitted.values),
       aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

ggplot(tibble("residuals" = ajuste$residuals,
                   "fitted" = ajuste$fitted.values^2),
       aes(x = fitted, y = residuals)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

De hecho fíjate que esa distribución leptocúrtica de nuestros errores hace que nuestra densidad de \(y\) sea un poco más apuntada que las simulaciones (que lo que hacen es usar el ajuste y añadirle el error simulando que fuese normal)

check_predictions(ajuste)
Warning: Maximum value of original data is not included in the
  replicated data.
  Model may not capture the variation of the data.

Fíjate como cambia al tratar mínimamente outliers

check_predictions(ajuste_outliers)

Ejercicio 6 (1/10)

¿Se cumplen las hipótesis? Argumenta porqué, más allá de lo que valga luego la bondad de ajuste, es importante que se cumplan.

dsadsa

Respecto a por qué es importante que se cumplan las hipótesis, podemos dar muchos motivos pero principalmente dos vistos en clase con ejemplos simulados:

  • No cumplir hipótesis implica no poder hacer inferencia: la recta solo nos sirve de manera muestral, para esta no muestra, no pudiendo establecer conclusiones del ajuste a la población. Esto, entre otras cosas, hace que no podamos interpretar ninguno de los p-valores y contrastes, ni que el modelo sirva más allá de esta muestra.

  • No cumplir hipótesis puede implicar que aunque la bondad de ajuste parezca buena, las predicciones sean bastante malas (por ejemplo, por heterocedasticidad)

Ejercicio 7 (1.4/10)

Evalua el modelo realizando un ANOVA, detallando y explicando todo lo que puedas cada uno de los elementos. Analiza la bondad de ajuste. ¿Qué % de información explica nuestro modelo? Detalla que implica respecto a la calidad del ajuste (en caso de que lo haga).

Para la significación gloval del modelo y el análisis de la varianza usamos anova()

anova_ajuste <- ajuste |> anova()
anova_ajuste
Analysis of Variance Table

Response: deathRate
                 Df  Sum Sq Mean Sq F value    Pr(>F)    
povertyPercent    1  343293  343293  540.47 < 2.2e-16 ***
Residuals      2411 1531400     635                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De aquí sacamos

  • que la suma de desviaciones al cuadrado de la variable estimada \(\hat{y}\) respecto a su media es de 343293, lo que hemos llamado \(SSR = 343293\)
  • que la suma de los residuos al cuadrado es \(SSE = 1531400\)
  • que \(SSE / (n-p-1)\) es aprox \(635\), que es lo que hemos llamado como \(\hat{\sigma}_{\varepsilon}^{2}\) (el estimador insesgado de la varianza residual), lo cual cuadra con la anterior salida donde Residual std error era de 25.2, que al cuadrado da 635.04
  • además en el contraste global \(H_0:\beta_1=\ldots=\beta_p=0\) (en este caso…\(p=1\)), obtenemos un p-valor muy pequeño –> rechazamos la hipótesis nula, el modelo tiene alguna predictora con un efecto lineal significativo respecot a la objetivo (en este caso, la única que tenemos)

Para el \(R^2\) podemos hacerlo de tres formas:

  • con summary(), que vemos que \(R^2 = 0.1831\)
ajuste |> summary()

Call:
lm(formula = deathRate ~ povertyPercent, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-121.490  -14.128    1.459   15.520  170.393 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    147.53777    1.44389  102.18   <2e-16 ***
povertyPercent   1.84645    0.07942   23.25   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 25.2 on 2411 degrees of freedom
Multiple R-squared:  0.1831,    Adjusted R-squared:  0.1828 
F-statistic: 540.5 on 1 and 2411 DF,  p-value: < 2.2e-16
  • con compare_performance() aunque solo tengamos un modelo
compare_performance(ajuste, ajuste)
Some of the nested models seem to be identical
# Comparison of Model Performance Indices

Name   | Model |   AIC (weights) |  AICc (weights) |   BIC (weights) |    R2 | R2 (adj.) |   RMSE |  Sigma
----------------------------------------------------------------------------------------------------------
ajuste |    lm | 22425.0 (>.999) | 22425.1 (>.999) | 22442.4 (>.999) | 0.183 |     0.183 | 25.192 | 25.203
  • de manera manual con el anova(), dividiendo \(SSR\) entre \(SST = SSR + SSE\)
anova_ajuste$`Sum Sq`[1] / 
  (anova_ajuste$`Sum Sq`[1] + anova_ajuste$`Sum Sq`[2])
[1] 0.1831195

El \(R^2\) nos dice que nuestro modelo solo consigue explicar el 18% aproximadamente, lo cual puede parecer poco si queremos predecir, pero a nivel médico, implica que solo con tener la variable pobreza de una persona podemos explicar el 18% de su mortalidad de cáncer a futuro (solo con una variable, en un tema tan complejo como el cáncer, not bad)

Ejercicio 8 (0.8/10)

Por último, utiliza el dataset test para, aplicando el modelo entrenado, predecir las observaciones de test. Con las observaciones de test construye un nuevo dataset de tres columnas: la variable objetivo, la variable predictora y la predicción. Haz lo mismo con train. Junta ambos datasets con una cuarta variable que nos diga si la observación era de train o de test.

Lo que hacemos primero es construir un dataset con la variable objetivo y predictora en test, y las predicciones usando predict() (ahora el newdata es literal el conjunto test que hemos apartado).

pred_test <-
  tibble("x" = test$povertyPercent,
         "y" = test$deathRate,
         "y_est" = predict(ajuste, newdata = test),
         "split" = "test")

Hacemos lo mismo con train y los juntamos (fíjate que he añadido esa cuarta variable para distinguir si las filas vienen de test o de train)

pred_train <-
  tibble("x" = train$povertyPercent,
         "y" = train$deathRate,
         "y_est" = ajuste$fitted.values,
         "split" = "train")
pred_data <- 
  bind_rows(pred_train, pred_test)
pred_data
# A tibble: 3,017 × 4
       x     y y_est split
   <dbl> <dbl> <dbl> <chr>
 1  13.5  194.  172. train
 2  16.1  198.  177. train
 3  41.9  248.  225. train
 4  20.6  220.  186. train
 5  22.1  214.  188. train
 6  20    174   184. train
 7  24.9  191.  194. train
 8  15.7  151.  177. train
 9  29.5  186.  202. train
10  28.2  216.  200. train
# ℹ 3,007 more rows

Ejercicio 9 (1/10)

Con el dataset anterior, visualiza de alguna manera la calidad de las predicciones en train y en test (idea: reales vs predicciones) en dos gráficos en el mismo ggplot (pista: la cuarta columna que has añadido antes la puedes usar en un facet). ¿En cuál ha acertado más el modelo? Calcula el \(R^2\) en el dataset de test (piensa como haciendo uso de los errores) y compáralo con el obtenido en la salida del lm() en train.

Respecto al R2 cuadrado podemos compararlos a mano sabiendo que es siempre ratio de varianza explicada.

R2_data <-
  pred_data |> 
  summarise(R2 = var(y_est)/var(y), .by = "split")
R2_data
# A tibble: 2 × 2
  split    R2
  <chr> <dbl>
1 train 0.183
2 test  0.179

IMPORTANTE: no podemos hacer un segundo ajuste en test, y calcular el R2. Lo importante de test es que vamos a evaluar la predicción con el ajuste ya estimado de train, ¡no uno nuevo! queremos “simular” que pasaría si ese modelo ya construido y entrenado, se lo pasamos a unos datos que el modelo nunca conoció

De ahí que salga algo peor en test que en train (lógicamente).

Para el gráfico, vamos a visualizar predicciones vs valores reales, en cada uno

ggplot(pred_data, aes(x = y, y = y_est)) +
  geom_point(aes(color = split), size = 2, alpha = 0.7) +
  geom_smooth(method = "lm") +
  facet_wrap(~split) +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

En ambos casos obtenemos algo que sabíamos: como modelo predictivo, es un modelo de poca capacidad predictiva.

¿Y si metiésemos más variables? Continuará…